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Abstract 



We present a dissipative algorithm for solving nonlinear wave-like equations 
when the initial data is specified on characteristic surfaces. The dissipative 
properties built in this algorithm make it particularly useful when studying 
the highly nonlinear regime where previous methods have failed to give a 
stable evolution in three dimensions. The algorithm presented in this work is 
directly applicable to hyperbolic systems proper of Electromagnetism, Yang- 
Mills and General Relativity theories. We carry out an analysis of the stability 
of the algorithm and test its properties with linear waves propagating on a 
Minkowski background and the scattering off a Scwharszchild black hole in 
General Relativity. 

I. INTRODUCTION 

When modeling nonlinear problems, dissipative algorithms have provided researchers 
with an extremely valuable tool since usually most non dissipative schemes fail to produce 
a stable evolution. More precisely, when using neutrally stable algorithms, instabilities 
often arise which spoil the evolution. The addition of artificial dissipation removes these 
instabilities by "damping" the growing modes of the solution, in a somewhat controlled way. 
Therefore, its inclusion in a discretization scheme provides a practical and economic way of 
achieving longer evolutions. 

The most widely used algorithms with this property are the family of Lax schemes ; 
whereby adding to the equation u t t = —au,x a term proportional to u >xx one obtains a stable 
discretization of the system that would otherwise be unstable. However, one might correctly 
ask whether this is not tantamount to solving a completely different problem. The beauty of 
these methods is that the proportionality factor depends on the discretization size, and in the 
continuous limit the approximation to modified PDE results in a consistent approximation 
to the original one. 
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Although there is much experience with these kind of schemes, most of the standard 
dissipative algorithms have been tailored for Cauchy initial value problems, where initial 
data is provided at one instant of time and evolved to the next instant by means of the 
evolution equation. However, in radiative problems, it is sometimes more convenient to 
choose a sequence of hypersurfaces adapted to the propagation of the waves, and therefore 
one adopts a foliation adapted to the characteristics of the PDE under study. 

In the present work, we present a new algorithm adequate for hyperbolic systems. 
The underlying strategy of the proposed algorithm is quite different from the conven- 
tional Cauchy-type methods. Rather, it is inspired by analytic concepts developed in the 
1960s 0-§J for studies of gravitational radiation in General Relativity and in their subse- 
quent numerical integrations^. The main features of this approach are the use of character- 
istic surfaces (for the foliation of the spacetime) and compactification methods (that enable 
the inclusion of infinity in the numerical grid) to describe radiation properties. Although 
evolution algorithms (for systems possessing some kind of symmetry) developed within this 
approach proved to be successful in the linear and mildly nonlinear regime |3|-§], they pro- 
duce unstable evolutions when applied to the general case; which shows the need for devising 
algorithms that could be applied in this situation. In the present work we present a new 
algorithm having dissipative properties, making it a valuable tool to study systems where 
strong fields might be present. 

In section 2 the details of the algorithm for the wave equation are presented and its 
stability properties discussed. Section 3 is devoted to treat a model problem which shows 
how the dissipative algorithm might be a useful tool for numerical investigations in General 
Relativity. Finally, in section 4 we describe two particular applications of this algorithm in 
the numerical implementation of General Relativity. 



For a detailed account of the developments in the characteristic formulation see || . 
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II. THE ALGORITHM 



Waves of amplitude g traveling in one spatial direction with unit speed obey the familiar 
equation 

g,tt-g,xx = o, (l) 

Which can be solved in the region R = {(t,x) ft > to,x £ R}, assuming g(t = t ,x) and 
9,t(t = to,x) are given. If, instead, one is interested in the solving the problem restricted to 
the region x £ [a, oo), boundary data must also be provided corresponding to g(t,x = a). 
The analysis of this problem can be described in a simple way in terms of the characteristics 
of this equation, which are given by (x — x ) = ±t through each spatial point x D . 

In particular, when solving Eq. (H) in the region IZc- The domain of dependence T>c of 
a point is given by V c = Sc f]T^c, with Sc naturally defined by the characteristics 

passing through (ti,Xi) as 

S c = {(t, x) such that t < t x and (t - ti) 2 - (x - x x ) 2 > 0} ; (2) 

and IZc is the region to the future of 

• the line t = to, 

or 

• the region defined by [a, oo) or x 6 [a, b] (where a £ R); in these cases, boundary 
conditions must be imposed at x = a (and x = b in the latter case). 

Suppose one introduces a coordinate system adapted to the characteristics by, say (u = 
t — x, r = x); then, Eq. (0) reduces to 

2£V - 9,rr = 0. (3) 

Further, one can then choose to foliate the spacetime by a sequence of characteristics ob- 
tained by holding the (retarded) time u = const. One can then define a characteristic initial 
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value problem, where Eq. (||) is solved provided that g{u = u a ,r) is given. (Note that 
9,u{u — u , r) need not be provided as in the Cauchy initial value problem). 



F{u) + G(u + 2r) (where F and G are smooth functions). Physically, F(u) describes waves 
propagating in the +r direction (outgoing waves) and G(u + 2r) describes waves propagating 
in the — r direction (incoming waves). Then, if one imposes the condition of pure outgoing 
waves, the solution must be of the form g = F(u); hence, along each characteristic the 
value of the function is constant. Notably, boundary data at r = can be given arbitrarily 
since purely outgoing waves at u = uq will not reach r = 0. More generally, boundary data 
consistent with the incoming waves must be prescribed at r = 0. 

It is important to note the domain of dependence for this problem. When solving Eq. @ 
in the region 7Z C , the domain of dependence (V c ) of a point (ui, r x ) is defined by V c = S c f] 1Z C 
where 



However, if the region 1Z C is chosen to be the future of the line u = Uq, V c extends indefinitely 
to the past. Therefore, the characteristic approach requires TZ C to have a boundary. Thus, 
one defines 1Z C as the region {u > u Q , r e [a, oo)) (with a > 0). Figure (|l|) illustrates the 
domains of dependence corresponding to each formulation. 

For hyperbolic systems with two or more spatial dimensions, the manner in which the 
characteristics determine the domain of dependence and lead to evolution equations is qual- 
itatively the same. Also, the use of coordinates adapted to them provide a tidy way for 
studying the system. For instance, in 3 dimensions, the wave equation is given by 



It is straightforward to check 




<S C = {(u,x) such that u <Ui and {u — Ui) 2 + 2{u — Ui){r — r x ) > 0}. 



(4) 




0. 



(5) 



which, in term of spherical polar coordinates (t, r, 8, 0) has the form 




(6) 



where L 2 denotes the angular momentum operator 
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sin(0) sin 2 (#)' {) 

Introducing coordinates (u = t — r,r, 8, 0), which defines a natural inner boundary at r = 0, 
Eq. @ takes the form 

2(r^), ur -(r^) irr = — . (8) 

r 

Thus, by defining g = r^fr and considering L 2 \l//r as a source term, Eq. (H) formally looks 
like the 1 dimensional system. Therefore, from now on we restrict our analysis to this latter 
case and extend our results to the 3 dimensional case in section 4. 

The formal integration of (|3]) proceeds by an integration in the r direction on each 
u = const surface and then evolve to the next level. This reformulates the integration in 
the characteristic formulation as an "evolution" in the radial direction and then another 
in the u direction (as opposed to the evolution of a "whole" instant of time to the next 
one typical of the Cauchy evolution). Hence, standard dissipative schemes intended for 
Cauchy-type evolutions (like the family of Lax algorithms) are not directly applicable in the 
characteristic formulation of the PDE and the addition of artificial viscosity to the system 
must be reformulated. 

In the numerical implementation of Eq. (|3|) a useful discretization was introduced in || . 
This scheme is basically a second order approximation based on finite difference techniques. 
Assuming the grid discretization is given by u n = nAu and = iAr, the derivatives may 
be discretized in the following way: 

n+1 n n+l n n _i_ n n 
in+1/2 _ 9j — 9j-l — 9j + 9j-l /q\ 

Mi-1/2 - AuAr > ^ 

n+1 O^.n+1 , „n+l , n o n n _j_ r .n 

,n+i/2 _ 9j - j9i-i + 9j-2 + 9i+i - Z 9i + 9i-l (m , 

9rr\i_ 1/2 ~ ■ l iU J 

The resulting scheme (which we will refer to as GIW) is a second order in time, second 
order in space accurate algorithm. Notably, the Von Neuman analysis shows that the GIW 
scheme has unitary amplification factor (i.e. a neutrally stable algorithm) independent of 
the values of Am and Ar. This would imply that the algorithm is unconditionally stable 



which is at first sight puzzling. This might be explained by the implicit local structure of the 
algorithm (since it involves 3 points at the upper time level) and, as such, a local stability 
analysis need not give a condition on the discretization size. Nevertheless, the algorithm is 
globally explicit as the evolution proceeds by an outward march from the origin. Hence, the 
algorithm does require the enforcement of the CFL condition to ensure that the numerical 
and analytical domains of dependence are consistent. 

The CFL condition for the system can be easily obtained. The field at grid point at 
(ui, ri), depends critically on the field value at {u\ — Aw, t\ + Ar) (since all the points where 
< r < Ti are trivially included in the discretization). The requirement for the numerical 
domain of dependence to include the analytical domain of dependence is Aw 2 — 2AuAr < 0; 
therefore, the CFL condition will be satisfied if Am < 2Ar. 

The GIW algorithm has been employed successfully in the characteristic formulation 
of General Relativity (G.R.) assuming either spherical symmetry [§,0; axisymmetry || 



or very small departures from spherical symmetry [[7|. However, when considering more 
general problems, as it is often the case with neutrally stable schemes, round off error or 
parasitic modes are enough to cause ripples in the solution which often lead to an unstable 
evolution. As stated earlier, adding dissipation to the PDE constitutes a way to alleviate 
this problem 0. We now show that a rather simple modification of (f|) can be used to obtain 
a consistent discretization with dissipative properties. 
We start by considering the following equation 

Ar 2 

2g,ur - 9,rr + 4 / 3e ^^^ r = ' 

(the 4/3 factor is included for convenience). A straightforward discretization of (|ll|) is 
obtained by the described approximation for g^ ur (^) and g >rr fllOD and by approximating the 
third derivative at the point (n, i — 1/2) as 

9,rrr\l 1/2 = ^jGfcl " 3tf + 3^ " 9?- 2 ) ■ (12) 

In analogy to the Lax method, the inclusion of this extra term leads to a consistent 
difference approximation of Eq. @; that is, the difference approximation converges formally 
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to the differential equation in the limit (Au, Ar) — > 0. In fact, it is straightforward to check 
that the resulting approximation is accurate of order {0(Ar 2 ),0(eAt)}^. An important 
feature of the resulting algorithm (which we shall call DA) is its dissipative features, which 
make it particularly useful. The stability properties of this algorithm can be easily obtained 
by introducing Fourier modes such that g = e SM e lfc i/ Ar . After some algebra one obtains 

S(i + 2asm(kAr/2)e- lkAr/2 ^j = z((l - e) + ^e(4 cos 2 (A;Ar/2) - 1)) - 2asm(kAr/2) e - ikAr/2 , 

(13) 

where S = e su and a = Aw/(4Ar). Therefore, the equation governing the growth of the 
solution's modes is 

' S '^ 1 + 3(l-^aW(L/2)) (- 2 + *W><*« + <") 

Now, since 4a(l — a) sin 2 (fcAr/2) < 1 (for a < 1/2) the scheme will be stable if, < e < 
3/2(1 — 2a). Moreover, \S\ < 1 as k — > n/Ar, indicating that the high frequency modes are 
effectively "damped", while \S\ — > 1 as k — > 0. 

The obtained discretization can also be thought of as an approximation to the original 
equation (|3]) (i.e. without the addition of the extra term) where the finite differencing of 
g^ ur includes 4 points on the nth level as 

in+l in 
I n+1/2 _ ^Vu-1/2 ~ P,r|j_i/ 2 

^Hi-i/2 - Au 

^ g- +1 - gg _ (i-^)(^-^-i) + ^r+i-e 2 )/3 (15) 

A-uAr AwAr ' 

which can be regarded as a weighted average of the derivatives at (n, i — 1/2) obtained from 
field values at the points {(n, i), (n, i — 1)} and {(n, % + 1), (n, i — 2)}. In the next section, 
we illustrate how this algorithm produces a stable discretization when the original strategy 
(corresponding to e = 0) fails. 



2 Contrary to the Lax-method which exhibits strict second order convergence in space and time. 
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III. APPLICATION IN A "TOY PROBLEM 



In this section we study the stability properties of an equation bearing close resemblance 
to the nonlinear evolution equation encountered in the characteristic formulation of General 
Relativity (which will presented in the section IV) 

2 G Mr — G rr — G G U G r . (16) 

In order to keep track of the nonlinearity of the equation, we introduce the parameter A 
(with A < 1), such that G = Xg; hence, Eq. ([16]) becomes 

2g, ur - g, rr = \ 2 g g tU g <r . (17) 



In particular, note that the principal part of Eq. flTTD corresponds to the wave equation. 
Also, it reduces in the the linear case (A 2 = 0), to the wave equation. Consequently; one 
might expect the GIW discretization to lead a stable scheme. 

However this is not the case, as can be demonstrated by the following analysis. First, in 
order to simplify the study of the stability of this nonlinear problem, we linearize the PDE 
with respect to the previous time step fl[ to obtain a more manageable equation. In this 
linearization, we approximate the values of g and g T with respect to the nth level, but g iU is 
centered in between the levels. The resulting finite difference approximation is 

A?/ 

+ ^ Nr +1 + 2 gtl - 9tl - 9? +2 + 2 g? - 9tx) = 

A 2 1(9? + 9tiM - 9tiM +1 + 9tl - 9? ~ 9ti) ■ (18) 

Finally, we introduce the Fourier modes g = e su e lk ^ Ar and solve for \S\ 2 , obtaining 

|5|2 i | 16 a A 2 sin(ir) 2 (l + cos(K)) ^ (ig) 

where a = Aw/(4Ar), K = kAr and 

D = 16 + A 4 (l + cos(fO) 2 - 8A 2 (asin 2 (fO + cos(i^)(l + cos(it))) 

+32a(l -a)(cos(AT) - 1) . (20) 
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It is not difficult to check that D is a positive quantity for 

o< Q .<I_^ + ^I±ZB^ (21) 

~ 2 8 8 V ; 

(which will be the case if the CFL condition is satisfied). Thus, the value of \S\ is always 
larger than 1, indicating that this discretization is unconditionally unstable! It is remark- 
able that the simple addition of some nonlinear terms, even when they do not change the 
equation's principal part, completely break an algorithm that would otherwise be stable. 

We now modify the way g jUr is discretized by introducing dissipation as dictated by the 
D.A. scheme, i.e. 

9,ur = (<?? - 97-x ~ ((1 - e)ff? - 9ti + <9? + x ~ 0?- 2 )/3)) ■ (22) 

With this simple modification, the value of \S\ 2 now becomes 

|g|2 = 1 + 4(00^-1)^ (23) 

where 

N = -4aA 2 (l + cos(K)) 2 

+e (A + e(cos{K) - 1) + A 2 cos(iT) (cos(iT) + 1) + 4a(cos(fT) - 1)) . (24) 

Since D > and cos(K) — 1 < 0, the condition for stability is that N > 0. Given a, this a 
condition on e, or vice versa. For instance, if a = 1/8 and A = 2~ 1 / 2 , the discretization will 
be stable if 0.3 < e < 1. On the other hand, if we choose A = 1 and e=l/4, iV>0 will be 
satisfied if a < |. Figure |3| shows the value of \S\ for different choices of e for a given A; the 
effect of the added dissipation can be clearly seen. 

IV. A PRACTICAL APPLICATION. THE CHARACTERISTIC FORMULATION 

OF G.R. 

When solving Einstein equations, one can take advantage of the coordinate invariance 
of the theory to simplify the modeling of a specific problem. In particular, one is free to 
choose a foliation of the spacetime better suited to the problem. 
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In the numerical implementation of G.R., the most common approach is to choose a 
foliation by spacelike hypersurfaces at constant times. In this approach, Einstein equations 
form a second oder PDE system for the intrinsic geometry of each surface and its embedding 
in the spacetime, the extrinsic curvature. Einstein equations are split into two distinct set of 
equations. One set consists of constraint equations that limit the possible configurations of 
the field variables on each hypersurface. The second set constitutes the evolution equations 
that determine the development of the field variables onto the next hyper surf ace^. 

The main drawback of the numerical implementation of the Cauchy formulation is the 
impossibility of having an infinite grid to completely cover the spacelike hypersurfaces. Thus, 
in practice, one chooses an exterior boundary in order to deal with a finite domain. This 
introduces further problems since special conditions on the boundary must be imposed in 
order to avoid reflections which spoil long term evolutions. Although in the 1 dimensional 
case there are sound method to achieve this (e.g. the Sommerfeld condition), in the general 
case, any local boundary condition still introduces reflections turning the task of obtaining 
long accurate evolutions into an almost impossible one. A related problem arising from an 
outer boundary at a finite distance is that the radiation can not be rigorously calculated. 

When studying gravitational radiation a more natural choice adapted to the wave prop- 
agation is to adopt a sequence of characteristic hypersurfaces to cover the spacetime. This 
approach is known as the characteristic formulation of G.R., pioneered by Bondi and 
Sachs The main ingredients of this formulation are the foliation of the spacetime 

by a sequence of characteristic hypersurfaces and the use of compactification techniques 
(which enable the inclusion of infinity in a finite grid) to rigorously describe asymptotic 
properties of radiation [||. The equations naturally split into hypersurface equations and 
evolution equations. We now outline the main aspects of the numerical implementation of 
this formulation (based on [j7|,|I2|; where a detailed description of the problem has been pre- 



3 For a complete description of this formalism see [11] 
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sented) and employ the constructed algorithm to discretize the PDE equations governing 
the evolution of the fields. 

A coordinate system is introduced by labeling the outgoing lightlike hypersurfaces with 
a parameter u; ach null ray on a specific hypersurface labeled with x A (A = 2, 3) and choose 
r as a surface area coordinate (i.e. surfaces at r = const have area 4nr 2 ). In the resulting 
x a = (u, r, x A ) coordinates, the metric takes the Bondi-Sachs form ||,[3| 

ds 2 = - {c 2p V/r - r 2 h AB U A U B ) du 2 - 2e 2l3 dudr - 2r 2 h AB U B dudx A + r 2 h AB dx A dx B , (25) 

where h AB h B c = $c anc ^ det{h AB ) = det(q AB ), with q AB a unit sphere metric. 
The metric components are re-expressed as 

h 22 = ^(^[J]+K), 
4 

h 2 3 = h 32 = -p2$*[J\, 

h S3 = ^(K-nJ}), 
U 2 = ^U[U], 

U s = ~m\ (26) 

where P = sec 2 (8/2) in standard angular spherical coordinates (9,(f)). Here, the metric is 
expressed in terms of two real {(3 and V) and two complex (U and J) variables (where 



K = vT-j- J J). The complex field J measures the departure of spherical symmetry of the 
surfaces given by r = const, and u = const; V represents the mass distribution of the 
system; j3 measures the expansion of the light rays and U measures the shift in the angular 
directions from one hypersurface to another (at constant r). 
The hypersurface equations are expressed as: 

Ar = W (27) 

U r = Fu[(3,J] (28) 
(r 2 Q), r = F Q {U,P,J] (29) 
V r = F v [Q,U,(3,J], (30) 
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where Q = r 2 e _2/3 ( J77 r + KU^ r ) which is introduced to deal with a first order system of 
hypersurface equations. The functions Fp, Fu, Fq and Fy involve derivatives taken only on a 
particular hypersurface M . Then, they can be easily integrated if J is known on M (assuming 
consistent boundary conditions are provided) in the following way. The integration strategy 
proceeds by first obtaining (3 from Eq. (P7|), then U from Eq. followed by the calculation 
of Q using Eq. ( p9|) and finally V using Eq. (|30|) . The evolution to the next hypersurface is 
prescribed by a first order (in time) equation for J that takes the form 

2(rJ) iUr -^(rJ) rT = rjfy{J r K-JK r ) + c.c.)+Fj[V,Q,U,(3,J] (31) 

where Fj involve derivatives on M only. 

A code that implements Einstein equations was written using (second order) finite differ- 
ence approximations. Angular and radial derivatives are approximated along the following 



lines HI2 



Angular derivatives. We follow the formalism given in [[U|,|l4]]. To expedite the nu- 
merical implementation of angular derivatives, instead of working with the standard 
spherical angular coordinates (9,<p), we work in stereographic coordinates: 

X A = (q iP ) = (tan(0/2)cos(0),tan(0/2)sin(0)) , (32) 
and angular derivatives are written in terms of the (complex differential) eth operators 



9 and 5 |15|,16|; for instance, 

dp _ 9/? + 9/3 
~dq ~ P 



(33) 



This allow us to employ a set of numerical techniques introduced in |L7[ which are 
specially tailored to: (i) handle the numerical approximation of angular derivative 
operators and (ii) deal with the fact that a single coordinate patch can not be used to 
smoothly cover the sphere. 

Radial derivatives. These are approximated via centered second order differences along 
each null ray (i.e. holding x A = const); for instance 
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/3r = A n _i + ArF /3 |r_ 1/2 . (34) 



The evolution equation deserves special consideration. Its discretization (in between 
levels n and n+1) is obtained using dissipation in the following way: schematically, it can 
be re-expressed as 

2g, ur - {V/r)g >rr = ^ {g >u g >r + c.c.) + Fj\J3, J, U, V}. , (35) 

where g = rJ. The function Fj can be straightforwardly approximated at each grid point 
(n + 1/2, i — 1/2) to second order accuracy. Then, in order to introduce dissipation in the 
algorithm, we proceed to consider a modified version (along the lines described in section 
II) 

g _ Ar 2 
2#,«r - (V/r)g >rr = {g, u g,r + c.c.) + Fj(p, J, U, V) + e-^-g >rrr . (36) 



We center the derivatives at the point (n+ 1/2, i — 1/2), as dictated by the DA scheme and 



obtain g, u \^i/2 by means of an iterative procedure. In the first iteration we set g^ u = g )U \i-i^ 2 



and get a first approximation of g? via the evolution equation. Then, we use this value to 
obtain a guess for g^^lj^ which is then used to get a better approximation for This 
procedure is repeated a sufficient number of times to ensure convergence. 

Unfortunately, when solving a 3 dimensional problem, the computational requirements 
of integrating from the origin (r = 0) are formidable. However, it is possible to start the 
integration from a finite value of r, assuming consistent values of /3, U, Q, V and J are 
known on this boundary (which is refered to as the worldtube boundary) as well as the 



value of J on an initial hypersurface [18 



To illustrate the usefulness of the presented algorithm, we apply it to model (i) the 
propagation of linear waves on a Minkowski background and (ii) the problem of scattering 
off a Schwarzschild black hole in 3 dimensions. 
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A. Linear Waves on a Minkowski background 



In the past, analytical solutions of linearized Einstein equations (in the characteristic 
formulation) have been found which describe waves propagating on a flat background 0. 
These solutions provide an important test bed for the algorithm, since the numerically 
obtained solutions must converge to the analytic values given by 

P = 0, V = r (37) 

with J and U obtained from a solution of the scalar wave equation by 



', = ™*. <*>) 



u, = -2*m+m. (39) 

In order to test the algorithm, we choose a solution of the wave equation in 3 dimensions, 
that represents an outgoing wave with angular momentum < I < 6 of the form 

rv 

<$> = {d z f — , (40) 

where d z is the z-translation operator. The resulting solution is well behaved above the 
singular light cone u — 0. 

Choosing initial data of very small amplitude (a ~ 10~ 9 ); we used these solutions to give 
data at it = 1 (with the inner boundary set at R — 1.5) and compared the numerical and 
exact solutions over time for different values of e. The computation was performed on grids 
of size N x equal to 41, 53, 65 (with the number of angular points N% = (N x — l)/2 + 5, and 
the ratio Aw/(4Ar) = 1/8). The L 2 norm of the error was calculated over the entire grid 
and plotted against time for different values of the dissipation parameter. Figure § shows 
the logarithm of the error in J vs. time (for runs with N x = 65). For e = the evolution 
is unstable, as can be seen by the exponential growth of the error. For e = 0.005 the 
instability appears at a later time, also with an exponential growth. However, for e = 0.05 
the run proceeds stably and the error remains under control. It is important to note that 

16 



the magnitude of the dissipation needed to achieve a stable run is very small and therefore 
the "damping" of the solution in not severe. 



B. Nonlinear scattering off a Schwarzschild black hole 

The characteristic initial value problem on an outgoing null hypersurface requires inner 
boundary conditions on the worldtube. Here we consider an example in which the inner 
boundary T consists of an ingoing nullcone (see figure |^). We adopt coordinates x A which 
follow the ingoing null geodesies and foliate V (chosen to correspond to ingoing r = 2m 
surface in a Schwarzschild spacetime) by slices separated by constant parameter u. In these 
coordinates, the Schwarzschild line element takes form 

(2m \ 
1 \du 2 - 2dudr + r 2 q AB dx A dx B . (41) 

The initial data correspond to setting J = as data on u = 0, with the boundary conditions 
(3 = U = Q = and V = r - 2m on T. 

We pose the nonlinear problem of gravitational wave scattering off a Schwarzschild black 
hole by retaining these boundary conditions on T but we choose null data at u = corre- 
sponding to an incoming pulse with compact support^, 



J(u = 0,r, x A ) 



(42) 

otherwise, 



where 2Yi,m is the spin-two spherical harmonic |T5|| . 

The code was run for different values of A under different choices of the dissipation 
parameter. In all cases, unstable evolutions resulted from the choice e = 0, however for 



4 Recall that the data on the initial hypersurface can be chosen freely in the characteristic formu- 
lation 101. 
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nonzero values of e the code ran without any stability problem as illustrated in figure |6| (for 
a run where A = 1, 1 — 2, m — 0). 

Yet, as expected of any dissipative algorithm, the solution decreases in amplitude with 
time. This highlights the need to carefully tune the value of e. Notwithstanding this fact, it 
is important to stress once again that this set of runs would not have been possible without 
dissipation. 

This problem was originally studied in the perturbative regime by Price ||19|| . There is 
no known analytic solution to the problem in the nonlinear regime and applying numerical 
methods is the only way to study it. The accuracy of the dissipative scheme can be assessed 
indirectly by inspection of the gravitational waves produced by the system. Gravitational 
waves can be described in terms of two polarization modes (refered to as plus and cross 
modes) |Tl| . However, when considering spacetimes with axisymmetry, the cross mode must 
vanish and this fact can be used to test the algorithm. Calculating the gravitational radia- 
tion is a rather involved problem that exceeds the scope of this work. A set of algorithms 
to numerically calculate the gravitational wave forms was constructed in the characteristic 
formulation in [T2| and tested under different situations. We used these algorithms in the 
present work to calculate the polarization modes for the choice e = Ar and with an axisym- 
metric pulse with / = 2,m = as the initial data. The cross polarization mode actually 
converges to zero in second order indicating an accurate discretization of Einstein equations, 
as can be seen in Figure [7|. 



V. CONCLUSION 

The algorithm described in this work represents a valuable tool for the study of nonlinear 
problems in the characteristic formulation. Its use enables long term evolution that would 
otherwise be impossible. Yet, there is still much room for improvement as the number of 
numerical techniques adapted to characteristic type evolutions is scarce (as opposed to the 
situtation in the Cauchy type evolution where one has at hand a great number of algorithms). 
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The variety of physical problems where propagating waves are to be described, highlights 
the need of further investigations on "characteristic" algorithms. 

Of particular interest is the application of these type of algorithms to the characteristic 
module constructed to model the collision of a binary black hole self-gravitating system. In 
this problem, it is imperative to have robust enough schemes capable of dealing with highly 
nonlinear fields. The complexity of the problem inspired the creation of the Binary Black 
Hole Grand Challenge Alliance, where a group of U.S. universities and outside collaborators 



are joining efforts to tackle the problem p0| . A strategy to study this problem is a "hy- 
brid" scheme that implements at the same time a Cauchy evolution (for the region near the 
black holes) and a Characteristic evolution (for the exterior region). This approach is called 
Cauchy- characteristic matching (CCM) pT| , p2| ,|8| , and in principle, its implementation man- 
ages to avoid the problems and to exploit the best features of each evolution scheme. CCM 
has been shown to work (and outperform traditional outer boundary conditions) in situa- 



tions where special symmetries were assumed [0,|T(| and its full 3 dimensional application 
in G.R. is currently under study. The characteristic code is one of the pieces of this bigger 
algorithm and the need for robust performance prompted this investigation. However, its 
use is not limited to G.R. Any hyperbolic system describing waves will have an evolution 
equation similar to Eq. @. The algorithm presented in this work should provide a useful 
tool in the numerical modeling of these systems. 
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FIGURES 



FIG. 1. Domains of dependence in the Cauchy and characteristic initial value problems. 

FIG. 2. The characteristic scheme for the exterior problem. Initial data is given on N and 
boundary data on T. 

FIG. 3. Plots of \S\ corresponding to different choices of e and A. A (A = 0, e = 0); B 
(A = 0.02, e = 0); C, (A = 0.02, e = 0.02) and D, (A = 0.02, e = 0.2) which illustrates how adding 
artificial dissipation ensures stability. However, as can be seen in D for a high value of e the 
damping of the high frequency modes might be severe. 

FIG. 4. The logarithm of \E\ = \ J num — J anal | (the numerical and analytic values of J) is shown 
for different values of e. For e = 0.05 the evolution is stable as opposed to the unstable evolutions 
correspondent to e = and 0.005. 

FIG. 5. Scattering off a Schwarzschild black hole. The bold dashed line illustrates the incoming 
pulse. 

FIG. 6. Plots of the field variable J at a representative angle vs. a compactified radial co- 
ordinate x = r/(l + r). The value of the mass is m = 0.5, the amplitude of the initial pulse is 
A = 5, R a = 4 and = 8. The runs correspond to different choices of e. The solid lines indicates 
the initial data at u = 1. Box (a) shows the run for e = 0, after a short time the obtained values 
diverge; (b) corresponds to the choice e = 0.005 showing a run that neither show signs of instability 
nor much damping of the pulse; (c) in turn corresponds to e = 0.02, although there is no sign of 
instability the solution has been damped considerably. 

FIG. 7. Convergence of the cross polarization mode to zero (in these runs e was chosen equal 
to Ar). The slope is 1.99, confirming second-order accuracy of the obtained wave form. 
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